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1.  Introduction 

We  discuss  the  queueing  model  GI  /M/c  In  which  customers  arrive  In  groups 
at  a c-server  unit  at  time  points  which  form  a renewal  process.  The  c ser- 
vers have  Independent  exponential  processing  times  with  the  same  rate  p. 

The  times  between  group  arrivals  have  the  same  distribution  F(*)  with  finite 
mean  a and  satisfying  F(CM-)  = 0.  The  successive  group  sizes  are  Independent, 
bounded  random  variables  with  probability  density  1 ^ v ^ K}  and  mean  n- 

Without  loss  of  generality,  we  require  that  6j^  > 0. 

It  Is  well-known  that  In  the  case  of  single  arrivals,  the  stationary 
queue  length  density  Immediately  prior  to  arrivals  exists  If  and  only  if 
c p a > 1 and  is  then  a modified  geometric  probability  density  in  which  the 
first  c terms  are  obtained  by  solving  a system  of  linear  equations.  We 
refer  to  the  book  by  D.  Gross  and  C.  M.  Harris  [6]  or  several  other  texts  on 
queueing  theory  for  a discussion  of  this  classical  result. 

The  extension  of  this  result  to  the  case  of  bounded  group  arrivals  appears 
to  be  new.  We  shall  obtain  the  stationery  queue  length  density  in  a modlf led 
matrix-geometric  form.  This  main  theorem  leads  to  a highly  efficient  algorithm 
and  yields  most  of  the  other  stationary  distributions  in  elegant  computable 
forms . 

The  queue  lengths  immediately  prior  to  arrivals  form  a Markov  chain  on 
the  nonnegative  integers.  The  transition  probability  matrix  P of  that  chain 
has  a form  which  depends  on  whether  K ^c  or  1 ^ K^c-1.  The  remainder 
of  this  section  Is  devoted  to  the  detailed  definition  of  the  matrix  P. 

We  define  the  quantities 
K 

(1)  ♦.(!)  - y (cpt)^'*'^,  for  j 1-K,  t ^ 0, 

^ v-max(l,-j) 

and 


r 


(2) 


a = 7 (Ji.(t)  dF(t),  for  j ^ -K. 

J Q J 


The  probability  generating  function  A*(z) 


K 


v=0 


a^_j^  z , is  given  by 


(3) 


A*(z)  = f*(cp-cviz)  ^ — — e z 
v=l  V 


K-v 


where  f*(*)  is  the  Laplace  - Stieltjes  tranform  of  F(’)»  The  mean  A*'(l-) 
is  given  by 


(4)  A*'  (1-)  = K - n + cpa. 

We  further  define  the  quantities  (t) , t ^0,  by 


(5)  it)  = (j) 


0. 

t 

f e 
0 


for  j r c, 
for  r < j ^ c. 


-cpt 


/ ^N^-c-l  ^Cv  -yj(t-T),,  -p(t-T),c-i, 

(cut)  cu  (.)  e Ml-e  ] •'dt, 

(r-c-1)!  ^ 

for  r > c,  0 < j c-1, 


and 


(6)  P 


ij 


v=l  0 


V i+v.j 


(t)  dF(t), 


for  i ^ 0,  0 ^ j c-1. 


The  transition  probability  matrix  P of  the  embedded  Markov  chain,  obtained 
by  considering  the  queue  length  immediately  prior  to  arrivals,  is  given  by 
(7)  For  K ^ c : 

P 
P 
P 
P 
P 
P 


00 

• • • Pn 

0,c-l 

d • • • 

-c 

Vk 

^-K 

0 

0 

0 

• • • 

10 

...  P-  , 

l,c-l 

^-c  ••• 

^2-K 

Vk 

®-K 

0 

0 

• • • 

a 

0 

• • • 

20 

. • • P , 

2, c-1 

^2-c  ••• 

^3-K 

®2-K 

®1-K 

"K 

30 

...  , 

3, c-1 

Ha  • • • 

3-c 

Vk 

®3-K 

Vk 

Vk 

«-K 

• * • 

40 

...  P , , 

4, c-1 

H < • • • 

4-c 

®5-K 

®4-K 

^3-K 

“2-K 

®1-K 

• • • 

50 

...  P,  , 

5, c-1 

H - • • • 

5-c 

®6-K 

^5-K 

®4-K 

®3-K 

^2-K 

• • • 

• 

• 

• 

• 

• 

• 

• 

• 

• 

I 

• 

• 

• 

9 

9 

• 

• 

• 

• 

• 

• 

m 
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For  1 < K < c-1: 


**00 

* * * 1 

0,c-l 

0 

0 

0 

• • • 

^c-K-1,0 

•••  ^c-K-l.c-l 

0 

0 

0 

• • 

I*  „ A 

P , , 

a 

0 

0 

* % * 

c-K,0 

...  c-K,c-l 

1 

^c-1,0 

• • • P - - 

^-1 

3^2  • • • 

3 

• • • 

^c,0 

. . . p , 

c ,c-l 

^0 

3 ^ • • • 

• • • 

p 

c+1,0 

"■  ^c+l,c-l 

ao  ... 

^2-k 

• • • 

^c+2,0 

•••  V2.C-1 

^2 

3^  • • • 

^3-k 

• • • 

The  first  c columns  are  given  by  Formula  (6) . The  remaining  columns  are  all 
obtained  by  successively  shifting  the  column  (a  , a,  • • • 1 downward  one 

place  at  a time. 

2.  The  Main  Theorem 
Theorem  1 

The  Markov  chain  P is  irreducible  and  aperiodic.  It  is  positive  recur- 
rent if  and  only  if  n < cya. 

For  K^c,  the  Invariant  probability  vector  3C  » (x^,  ...),  parti- 

It 

tloned  into  K -vectors,  is  given  by  Xj^  ■ JC^R  , for  k ^ 0.  The  matrix  R is 
strictly  positive  and  its  spectral  radius  is  less  than  one. 


Defining  the 

Kx 

K matrices 

A , n : 
n 

L 0,  by 

^-K 

0 

...  0 

^nK-K 

* • • ®nK-2K+l 

®1-K 

«-K 

...  0 

®nK-K+l 

• • • ®nK-2K+2 

0 “ 

• 

• 

• 

®-l 

• 

• 

• 

«-2 

...  a_^ 

1 

• 

• 

®nK.l 

• 

• 

• 

• • • ®nK-K 

n > 1 
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the  matrix  R is  the  minimal  solution  of  the  matrix  equation 
(9) 


R = 


n=0 


r"a  , 

n 


in  the  set  X of  nonnegative  matrices  of  order  K with  spectral  radius  less 
than  one.  The  matrix  R is  the  unique  solution  of  (9)  in  the  set  X. 

Defining  the  K x K matrices  B^,  n ^ 0,  by  partitioning  the  first  K 
columns  of  the  matrix  P for  K ^ c into  K x K blocks,  the  matrix 

(10) 


B[R]  = 2L.  ^\> 

n=0  " 


is  a strictly  positive,  stochastic  matrix.  The  vector  ^ is  the  unique 
vector  which  satisfies 


XQ(I-R)'^e  = 1, 


(11) 

where  je  =•  (1,1 1)  ’• 

For  1 ^ c-1,  we  partition  the  Invariant  vector  as  (3^,  X2»***)» 

where  Xj^,  k 5^  1,  has  K components  and  c-K  components.  The  vectors  jc^, 

k-1 

k > 1,  are  given  by  Xj^  * x^^  > where  the  matrix  R is  defined  in  the  same 

manner  as  above. 

We  partition  the  first  c-r  columns  of  P for  K ^ c-1,  into  a (c-K)  x 
(c-K)  block  Bqq  and  the  remainder  into  K x (c-K)  blocks  for  n ^1. 

The  next  K columns  are  partitioned  into  a first  (c-K)  x K block  and  the 


remainder  into  K x K blocks  for  n ^ 1- 


The  c X c matrix  B[R],  given  by 


(12) 


B[R] 


00 


n«l  n“l 


01 

r"-^b 


nl 


is  irreducible  stochastic.  The  c-vector  (Jt^.Xj^)  is  the  unique  vector 


satisfying 


B[R].  2So®  + (I-R)  e = 1. 


Proof 


The  matrix  equation  (9)  was  examined  in  detail  for  a more  general 
situation  in  [11].  We  shall  draw  on  results  proved  there,  without  repeating 
the  lengthy  proofs. 

00_ 

The  matrix  A = A^  is  readily  seen  to  be  a strictly  positive 

circulant  stochastic  matrix.  Its  invariant  probability  vector  ^ is  given  by 
1 00 

tr  = K e'.  The  vector  B = r nA  e satisfies  n B = K ^ A*' (1-)  = 

— — — n=l  n—  

K ^(cpo  + K-n) . This  equality  was  proved  in  [12],  Lemma  1. 

In  [11],  we  showed  that  a unique  matrix  R with  the  stated  properties 

exists  if  and  only  if  £ > 1,  or  equivalently  cpa  > q.  If  cpa  q,  we 

know  that  the  queue  cannot  be  stable.  Since  the  matrix  P is  clearly  irre- 
ducible and  aperiodic,  we  establish  positive  recurrence  by  proving  that  P 
has  a strictly  positive  invariant  probability  vector  if  cpa  > q. 

For  K ^ c,  the  remaining  statements  all  follow  by  application  of  Theorem 
2 in  [11].  Since  the  matrices  R and  B^,  n ^ 0,  are  positive,  so  is  B[R] 
and  hence  the  vector 

For  1 ^ c-1,  the  matrix  B[R]  is  positive,  except  for  a triangular 

corner  on  the  upper  right  hand  side.  Such  a matrix  is  clearly  irreducible. 


In  order  to  verify  that  B[R]  is  stochastic,  it  suffices  to  check  that 


n“l 


n-1 


6 


The  left  Invariant  vector  of  B[R],  normalized  so  that 

+ Xj^Cl-R)  = 1,  is  uniquely  determined  and  positive.  It  is  now  a routine 

2 

matter  to  verify  that  the  positive  probabi.’ ' ty  vector  x = (Xq  »2L]^  «•••) 

is  the  invariant  vector  of  P. 

Corollary  1 

For  K ^ c,  the  stationary  mean  queue  length  including  the  customers 
in  service,  prior  to  arrivals  is  given  by 


(lA) 


= y~  XqR'^  (vKe  + u)  = 2^0  (I-R)  ^ H + Kx^  (I-R)”^Re, 


where  u = (0,  1,  2,  K-1) ' . 

For  1 1.  K £ c-1,  the  corresponding  formula  is 


(15) 


00 

= XqH.'  + ^ ^ f (c-K)£  + (v-1) 


+ Xj^Cl-R)  ^ [u  + (c-K)^]  + Kx^  (I-R)  ^ R£, 

where  u'  = (0,  1,  c-K-1) ' . 

Corollary  2 

If  the  customers  in  arriving  groups  are  numbered  1,  2, 


, so  as  to 


specify  the  order  in  which  they  will  enter  service,  the  stationary  waiting 
time  distribution  W^(*),  1 ^ r ^ K,  of  the  r^^  customer  in  a group  is  given 
by 


U(x)  + ? 


(16)  W (x)  = X.  'U(x)  . . lu,. 

r l<c-r  i i’=max(0  ,c-r+l)  i 1+r-c 


X.  E, 


(cp,x). 


for  X ^ 0,  where  U(*)  is  the  degenerate  distribution  and  Ej^(cu,*)  is  the 
Erlang  distribution  of  order  k with  parameter  cp. 

Proof 

With  probability  x^,  the  r^^  customer  in  the  arriving  group  finds 
1+r-l  customers  ahead  of  him.  If  1+r-l  f.  c-1,  he  does  not  have  to  wait. 


7 


If  i+r-1  ^ c,  he  waits  until  the  number  of  customers  ahead  of  him  drops  to 

c-1,  or  equivalently  until  i+r-c  services  have  been  completed.  Since  during 

that  time  all  c servers  are  busy,  the  distribution  of  the  time  to  complete 

i+r-c  services  has  an  , (cu,*)  distribution. 

i+r-c 

Corollary  3 

The  stationary  distribution  W*(*)  of  the  server  backlog,  immediately  after 
the  arrival  of  a group  is  given  by 


(17) 


«*(X)  h)  + 

r=l  i<c-r 


1 j _ 


r=l  i=max(0,c-r+l) 


where  ^ j (^t)  “ (1-e  , for  x ^ 0. 


J 


Proof 


Conditioning  on  the  size  of  the  arriving  group,  we  see  that  if  r cus- 
tomers join  and  there  are  at  most  c-r  customers  present,  the  backlog  con- 
sists of  the  time  until  i+r  customers  are  processed.  The  distribution  of 
that  time  is 

If  r customers  Join  and  there  are  already  i present,  so  that 

i+r  > c,  the  backlog  consists  of  the  time  until  only  c customers  remain, 

plus  the  time  required  to  serve  those  c customers  by  c parallel  servers. 

This  leads  readily  to  the  second  term, 

3.  The  Queue  Length  in  Continuous  Time 

In  this  section,  we  assume  that  F(*)  is  not  a lattice  distribution. 
Analogous,  but  cumbersome  formulas  hold  in  the  lattice  case. 

By  we  denote  the  vector  of  probabilities  yj , j 0,  where  yj  is  the 
probability  that  j customers  are  present  at  time  t in  the  stationary  version 
of  the  queue.  We  shall  find  it  convenient  to  partition  the  vector  ^ in  the 


8 


L 


same  manner  as  the  vector  x was  partitioned  above.  We  shall  again  need  to 
distinguish  between  the  formulas  for  the  cases  K ^ c and  1 ^ K ^ c-1. 

It  is  also  convenient  to  introduce  a stochastic  matrix  P which  has  the 


same  general  form  as  the  matrix  P,  defined  in  Formula  (7),  except  that  the 
corresponding  entries  are  defined  by 

(18) 


K 


= / ‘t'4  (t)  ^ dt,  for  j 1 -K 

J Q J 


/ 


P,.  = ^0  1 . (t)  - — ^ dt,  for  i > 0,  0 < j < c-1. 

ij  v=l  V 0 l+v,j  a — — — 


We  shall  partition  the  matrix  P in  the  same  manner  as  in  the  proof  of 
Theorem  1 and  the  corresponding  matrices  in  the  partition  will  be  denoted  by 
the  same  letters,  but  distinguished  by  tildes. 

By  a classical  application  of  the  Key  renewal  theorem,  the  probabilities 


Yj t j ^ 0,  are  related  to  the  probabilities  x , j ^ 0,  by 
Jr™  j 


(19)  y 


j 


Tie  y " 

r=l  r i=max(0, j-r)  a 0 


J (u)  [l-F(u)  ]du,  j i 0, 


where  for  u > 0, 


(u)  = 0 


for  j > i ^ 0, 


e (cuu)^  for  i ^ j ^ c, 

(l-j)l 


(‘)  (1-.-““)^^  .-"J".  tor  0<  j < 1<  0. 


for  1 > c > j ^ 0. 

Formula  (19)  is  not  particularly  enlightening,  but  shows  upon  inspection 
that 

(20)  y - X P. 

By  considering  the  partitioned  forms  of  the  matrix  P,  we  obtain  computa- 
tionally useful  formulas. 


Theorem  2 


For  K ^ c,  we  obtain 

00 

(21)  Jo  = Xq  i[Rl  = Xq  ZI 


n»0 


^ ^ A[R],  for  ill. 


where  A[R)  = r”a  . 


n=0 


For  1 <■  K <■  c-1,  we  obtain 


(22)  ^ ®00  BnO 

n=l 


il  ■ ^0  ®01  ^ il  ®nl 

n=l 


^ = ^1 


.1-2 


X,  R A[R], 


for  1 > 2. 


Proof 


By  direct  computation  and  by  application  of  Theorem  1. 

Corollary  4 

For  K 1 c,  the  mean  queue  length  at  an  arbitrary  time  point  in  the 
stationary  version  of  the  queue  is  given  by 

(23)  4 “ *0  (I-R)"^A[R]^  “ (I-R)”^  A[R]e, 

where  u^  ■ (0,  1,  K-1) ' • 

The  corresponding  formula  for  the  case  1 1 K 5. 

(24)  L2  - 2So  BqO  ^ ®01^  ^®01^ 


-1 


+ (c-K)  Xj^  Bj^lRl®.  212  ii  2I2  (^“B)  A[R]  u, + 


Kxj^  (I-R)”^  A[R]e  + KXj^  (l-R)“^  AtR]e, 


f 


I 


10 


where  u = (0,1, . . . ,K-1)  ’ , u'  = (0,1, . . . ,c-JC-l) , B-(R]  - 5^  R"~^B 

u , nO 

«. 

and  B, [R]  . ^ ~ 

i , nl 

n=l 

4.  Numerical  Methods 

We  first  note  that  the  matrices  P and  P,  for  a general  value  of  K 
are  simply  related  to  the  corresponding  matrices  for  K “ 1.  Let  Q and  Q 
be  the  transition  probability  matrices  corresponding  to  the  case  K * 1, 
and  let  for  any  Infinite  stochastic  matrix  X with  Index  set  (0,1,2,...}  , 
the  "shifted"  matrix  AX  be  defined  by 

(25)  ^^^IJ  “ ^i+l,j*  forl^O,  j^O, 

then  It  is  easy  to  see  that 


(26) 


P “ e^Q  + 0,AQ  + ...  + e^A*^"^Q, 

.K-1; 


P - + e2AQ  + ...  + 0jjA  Q. 

It  Is  therefore  sufficient  to  discuss  the  computation  of  the  entries  of 
the  matrices  Q and  Q.  Except  for  the  first  c columns,  all  entries 
of  Q and  Q Involve  only  the  quantities 


(27)  b - / -iSHi 


dF(t) , dt. 


for  V > 0. 


A simple  partial  Integration  shows  that 

V 

(28)  b * (1  - y b , . for  v ^ 0, 


b - — (1  

V cpa  j 


so  that  only  the  quantitltes  {b^}  remain  to  be  found.  In  general,  t?-'®  latter 
require  a careful  numerical  Integration.  If  F(‘)  has  a density,  quadrature 
formulas  using  the  properties  of  Laguerre  polynomials  may  be  profitably 
Implemented.  If  F(*)  Is  a discrete  distribution,  the  probability  density 
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{b  } is  a mixture  of  Poisson  distributions  and  the  b can  be  evaluated  by 
a simple  recursive  scheme  using  the  recurrence  relation  for  the  Poisson 
probabilities.  In  the  important  case,  where  F(*)  is  a distribution  of  phase 
type  [8],  the  probabilities  can  also  be  evaluated  by  a numerically 

stable  recursive  method,  which  does  not  involve  any  numerical  integrations 
[9,10]. 

In  order  to  compute  the  first  c columns  of  Q,  we  first  note  that  the 
rows  labeled  0,...,c-l,  have  1, . . . ,c  positive  elements  respectively.  More- 
over the  row-sums  are  equal  to  one.  These  entries  may  be  computed  by  routine 
numerical  integrations. 

Next,  we  observe  that 


r-c+1 


= 1 - 


b = CPU  b . , , 
V r-c+1 


for  r > c. 


Also,  for  0 j c-1. 


r=c  j=0 


= cpa 


These  relations  can  be  used  as  accuracy  checks  and  also  to  define  truncation 
rules  in  procedures  to  evaluate  the  quantities  Qj.j»  for  r ^ c,  0 ^ J ^ c-1. 

In  general,  these  quantities  may  be  computed  by  solving  the  system  of 
differential  equations 

(31)  (t)  “ -VJ  ♦r+ij  (b)  + P(j+1)  ^r+l,j+l^’^^  ’ 0 1 J 1 c-1, 

numerically,  after  truncation  at  a sufflcently  high  index  r.  We  note  that 

•“CUt 

the  system  is  Inhomogeneous,  since  (t>  (t)  is  equal  to  e (cpt)  , 

(r-c+1)! 

for  r > c.  The  Initial  conditions  are  t .(0)  * 0,  for  r > c,  0 < j < c-1. 

— r+l,j  — — — 


The  probabilities  are  then  obtained  by  a progressive  numerical  integration 
method,  such  as  Simpson's  rule. 
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A routine  partial  Integration  yields  that 


(32)  . 

The  quantities 


may  be  found  by 


for  1 ^ j £ c-1,  r ^ c. 


(33) 


‘rO 


c-1  r-c 


for  r > c. 


The  preceding  discussion  applies  to  any  distribution  F(*)  and  lists 

computational  shortcuts,  which  are  useful  In  all  cases.  If  F(*)  Is  a 

probability  distribution  of  phase  type.  It  Is  possible  to  compute  all 

entries  of  Q and  Q by  recursive  schemes,  which  do  not  require  any  numer- 
ical Integrations  whatever.  The  details  of  this,  as  well  as  extensive 

numerical  examples  will  be  discussed  elsewhere. 

The  next  major  step  In  the  algorithm  Is  the  computation  of  the  matrix  R. 

The  equation  (9)  may  be  quite  efficiently  solved  by  successive  substitution 

In 


(34)  R = J_r"  a (I-A.)"^, 

nyi 

starting  with  R > 0.  Extensive  numerical  experience,  with  matrices  of  order 
as  high  as  fifty.  Indicates  that  convergence  Is  rapid,  except  for  queues 
which  are  close  to  critical. 

The  matrix  B[R]  may  be  routinely  computed.  Since  B[R]  is  stochastic, 

the  row-sums  of  the  computed  matrix  should  be  close  to  one.  Row-sums  ex- 
-9 

ceedlng  1-10  can  usually  be  attained  without  excessive  computational 
efforts. 

r 

A further  accuracy  check  la  obtained  by  verifying  that  the  matrix 
B[R]  + (I-R)  ^ A[R]  Is  stochastic.  In  practice,  a large  number  of  terms 
of  the  probability  densities  x and  as  well  as  the  distribution  of  waiting 


! 

times  and  a variety  of  moments  may  be  computed  at  a small  cost  and  to  an 
adequate  degree  of  accuracy. 

5.  Related  Queueing  Models 

1: 

The  literature  on  the  GI/M/c  queue  and  Its  various  generalizations  Is 

1 1 

fairly  extensive.  A useful  set  of  references  Is  given  In  the  paper  [2]  by  | 

U.N.  Bhat,  which  deals  with  transient  results  In  the  case  of  single  arrivals.  ' 

I 

F.  G.  Foster  [4]  and  F.  G.  Foster  and  A.  G.  A.  D.  Perera  [5]  deal  with 

' i 

group  arrivals  of  fixed  group  sizes  for  the  single  server  case.  It  may  be 

i! 

verified -that  the  roots  Inside  the  unit  disk  of  the  equation  (1)  In  [4]  are  ; 

k 'i 

the  eigenvalues  of  the  matrix  R,  corresponding  to  the  GI  /M/1  queue,  discussed  I 

there . 

The  extensions  to  bulk  service,  treated  by  P.  B.  M.  Roes  [13]  and  K.  ‘ 

J} 

Shyu  [14],  affect  only  the  analytic  expressions  for  the  entries  of  the  matrix  f 

t' 

P,  but  not  the  general  nature  of  the  partitioning  of  P.  The  latter  determines  \ 

the  form  of  the  steady-state  vector,  so  that  the  bulk  arrival-bulk  service 

extension  of  the  GI/M/c  queue  may  be  treated  by  routinely  adapting  the  method  ' 

1 

given  here . I 

I 

The  algorithmic  papers  by  F.  S.  Hllller  and  F.  D.  Lo  [7],  0.  S.  Yu 
[15],  and  D.  M.  Avis  [1]  deal  with  Erlang  service  times  and  heterogeneous 

i 

servers,  but  Involve  single  arrivals.  Although  also  In  this  case,  a matrix-  | 

geometric  characterization  of  the  steady-state  vector  may  be  given,  the  | 

dimension  of  the  matrix  R Is  so  large  that  It  does  not  provide  us  with  a 
feasible  computational  alternative  to  the  method  proposed  by  these  authors. 

Finally  we  note  that  the  GI  /M/c  queue  with  bounded  group  sizes  for  the 
arrivals,  may  also  be  considered  as  a particular  case  of  an  SM/M/c  queue  In 
which  some  of  the  sojourn  time  distributions  of  the  seml-Markovlan  arrival 
process  are  degenerate.  Details  of  this  construction  may  be  found  In  E. 
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^Inlar  [3,  p.  378].  The  matrix-geometric  characterization,  given  in 
Corollaries  3 and  4 of  Neuts  [11],  of  the  steady-state  queue  length  distri- 
bution for  the  SM/M/1  queue  may  be  generalized  in  a direct  manner  to  the 
SM/M/c  queue.  This  approach  would  lead  to  the  same  computational  results 
via  different  algorithmic  steps  and  there  may  be  some  merit  in  an  empirical 
comparison  between  the  two  methods . 
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